# Load Package ------------------------------------------------------------
library(xtable)
library(survey)

# Load Data ---------------------------------------------------------------
# Set file path
setwd("")

load("drop_in_ocean_data.rdata")


# Table -------------------------------------------------------------------
# Combined and formatted outside of R
# Age
age_out <- cbind(round(prop.table(table(cjs.df$ageDem[which(cjs.df$black == 1)]))*100),
                 round(prop.table(svytable(~ageDem, subset(d.bw.stk, black == 1)))*100),
                 round(prop.table(table(cjs.df$ageDem[which(cjs.df$black == 0)]))*100),
                 round(prop.table(svytable(~ageDem, subset(d.bw.stk, black == 0)))*100))

print(xtable(age_out, 
             caption = "Age Distribution",
             digits = 0,
             align = c("clccc")),
      include.rownames=F)

# Education
educ_out <- cbind(round(prop.table(table(cjs.df$edu5[which(cjs.df$black == 1)]))*100),
                  round(prop.table(svytable(~edu5, subset(d.bw.stk, black == 1)))*100),
                  round(prop.table(table(cjs.df$edu5[which(cjs.df$black == 0)]))*100),
                  round(prop.table(svytable(~edu5, subset(d.bw.stk, black == 0)))*100))

print(xtable(educ_out, 
             caption = "Education Distribution",
             digits = 0,
             align = c("clccc")),
      include.rownames=F)

# Female
fem_out <- cbind(mean(cjs.df$man[which(cjs.df$black == 1)])*100,
                 svymean(~man, subset(d.bw.stk, black == 1))*100,
                 mean(cjs.df$man[which(cjs.df$black == 0)])*100,
                 svymean(~man, subset(d.bw.stk, black == 0))*100)

print(xtable(fem_out, 
             caption = "Sex Distribution",
             digits = 0,
             align = c("clccc")),
      include.rownames=F)

# Income
inc_out <- cbind(round(prop.table(table(cjs.df$demInc[which(cjs.df$black == 1)]))*100),
                 round(prop.table(svytable(~demInc, subset(d.bw.stk, black == 1)))*100),
                 round(prop.table(table(cjs.df$demInc[which(cjs.df$black == 0)]))*100),
                 round(prop.table(svytable(~demInc, subset(d.bw.stk, black == 0)))*100))

print(xtable(inc_out, 
             caption = "Income Distribution",
             digits = 0,
             align = c("clccc")),
      include.rownames=F)